SUBROUTINE WRTCLIM
USE PARKIND1  ,ONLY : JPIM     ,JPRB,  JPRD
USE YOMHOOK   ,ONLY : DR_HOOK, JPHOOK  ,LHOOK
! (C) Copyright 2000- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!**** *WRTCLIM*  - writes fixed climate fields to NetCDF files

!     Purpose.
!     --------
!        writes out fixed climate fields

!**   Interface.
!     ----------
!        *CALL* *WRTCLIM

!     Explicit arguments :
!     --------------------


!        Implicit arguments :
!        --------------------

!     Method.
!     -------
!        See documentation

!     Externals.
!     ----------
!        None

!     Reference.
!     ----------

!     Author.
!     -------
!        Bart vd Hurk (KNMI)

!     Modifications.
!     --------------
!        original: 14/7/2000
!        I. Ayan-Miguez (BSC) June 2023: Added surface spatially distributed parameters
   
USE YOMLUN1S , ONLY : NPOSCLM  ,RMISS, NULOUT
USE YOMLOG1S , ONLY : LWRCLM   ,NACCUR   ,NDIMCDF
USE YOMGC1S  , ONLY : LMASK
USE YOEPHY   , ONLY : LEVGEN
USE YOESOIL1S, ONLY : RDAW     ,RWPWP    ,RWSAT , RWCAP
USE YOMDPHY  , ONLY : NLAT     ,NLON     ,NCSS,  NPOI, NPOIOFF &
     &               ,NLALO    ,NSOTY , YDSURF, NPOIP
USE YOMGPD1S , ONLY : VFCVL    ,VFCVH    ,VFTVL    ,VFTVH, VFCO2TYP, S2RVCOVH2D, S2RVCOVL2D,&
 & S3RWPWPM3D, S3RWSATM3D, S3RWCAPM3D, S3RVROOTSAH3D, S3RVROOTSAL3D
USE MPL_MODULE
IMPLICIT NONE

!* Netcdf Interface
INTEGER NVARID,NPOS,IERR
INTEGER ISTART3(3)
INTEGER ICOUNT3(3)
REAL*4,ALLOCATABLE :: ZOUTPUT3S(:,:)
REAL*8,ALLOCATABLE :: ZOUTPUT3(:,:)

!* Local variables
REAL(KIND=JPRB),ALLOCATABLE :: ZVALUE(:),ZBUF(:)
INTEGER(KIND=JPIM) IPOS,ITYPL,ITYPH
LOGICAL LPOS
INTEGER(KIND=JPIM) :: IX1,IX2,IY1,IY2,NLEN,IST,IEN,JK,JL,ISTP,IENP
REAL(KIND=JPRB) :: ZCVL,ZCVH
INTEGER(KIND=JPIM) :: MYPROC, NPROC
REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

#include "surf_inq.h"
#include "netcdf.inc"

IF (LHOOK) CALL DR_HOOK('WRTCLIM',0,ZHOOK_HANDLE)

MYPROC = MPL_MYRANK()
NPROC  = MPL_NPROC()

ISTP   = NPOIOFF(1)+1
IENP   = NPOIOFF(NPROC+1)
WRITE(NULOUT,*) 'WRTCLIM'


IF(NDIMCDF == 2)THEN
  IX1=1
  IX2=NLON
  IY1=1
  IY2=NLAT
  NLEN=(IX2-IX1+1)*(IY2-IY1+1)
  IST=IX1+(IY1-1)*NLON
  IEN=IST+NLEN-1
ELSE
  IST=1
  IEN=NLALO
  NLEN=IEN-IST+1
ENDIF
ALLOCATE (ZOUTPUT3(NLEN,NCSS))  !INFLATED OUTPUT ARRAY
ALLOCATE (ZOUTPUT3S(NLEN,NCSS)) !INFLATED OUTPUT ARRAY

ALLOCATE (ZBUF(NLEN))
ZBUF(:)=RMISS
!* ===============================================================
!     START OUTPUT
!* ===============================================================

IPOS=NPOSCLM
LPOS=LWRCLM

IF(NDIMCDF == 2)THEN
  ISTART3(1)=IX1
  ISTART3(2)=IY1
  ISTART3(3)=1
  ICOUNT3(1)=IX2-IX1+1
  ICOUNT3(2)=IY2-IY1+1
  ICOUNT3(3)=NCSS
ELSE
  ISTART3(1)=IST
  ISTART3(2)=1
  ICOUNT3(1)=NLEN
  ICOUNT3(2)=NCSS
ENDIF

ALLOCATE (ZVALUE(NPOI))

IF(LWRCLM)THEN

!* -- The fixed variables
!******************************************************************
  IF (LEVGEN) THEN
    CALL SURF_INQ(YDSURF,PRDAW=RDAW)

  ELSE
    CALL SURF_INQ(YDSURF,PRDAW=RDAW,PRWSAT=RWSAT,PRWCAP=RWCAP,PRWPWP=RWPWP)
  ENDIF

  NPOS = NPOSCLM

!     soil depth
  IF( MYPROC == 1 ) NVARID = NCVID(NPOS,'SoilThick',IERR)
  DO JK=1,NCSS
    ZVALUE(:)=RDAW(JK)
    CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZVALUE(:),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTCLIM:")
    IF( MYPROC == 1 ) THEN
      ZOUTPUT3(:,JK)= UNPACK(ZBUF,LMASK(ISTP:IENP),RMISS)
    ENDIF
  ENDDO
  IF( MYPROC == 1 ) THEN
    IF(NACCUR == 1)THEN
      ZOUTPUT3S(:,:)=ZOUTPUT3(:,:)
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT3S,IERR)
    ELSE
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT3,IERR)
    ENDIF
  ENDIF

!     field capacity
  IF( MYPROC == 1 ) NVARID = NCVID(NPOS,'SoilFC',IERR)
  DO JK=1,NCSS
    IF(LEVGEN)THEN
      DO JL=1,NPOI
        ZVALUE(JL)=S3RWCAPM3D(JL,JK)
      ENDDO
    ELSE
      ZVALUE(:)=RWCAP
    ENDIF
    CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZVALUE(:),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTCLIM:")
    IF( MYPROC == 1 ) THEN
      ZOUTPUT3(:,JK)= UNPACK(ZBUF,LMASK(ISTP:IENP),RMISS)
    ENDIF
  ENDDO
  IF( MYPROC == 1) THEN
    IF(NACCUR == 1)THEN
      ZOUTPUT3S(:,:)=ZOUTPUT3(:,:)
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT3S,IERR)
    ELSE
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT3,IERR)
    ENDIF
  ENDIF

!     wilting point
  IF( MYPROC == 1 ) NVARID = NCVID(NPOS,'SoilWilt',IERR)
  DO JK=1,NCSS
    IF(LEVGEN)THEN
      DO JL=1,NPOI
        ZVALUE(JL)=S3RWPWPM3D(JL,JK)
      ENDDO
    ELSE
      ZVALUE(:)=RWPWP
    ENDIF
    CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZVALUE(:),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTCLIM:")
    IF( MYPROC == 1 ) THEN
      ZOUTPUT3(:,JK)= UNPACK(ZBUF,LMASK(ISTP:IENP),RMISS)
    ENDIF
  ENDDO
  IF( MYPROC == 1) THEN
    IF(NACCUR == 1)THEN
      ZOUTPUT3S(:,:)=ZOUTPUT3(:,:)
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT3S,IERR)
    ELSE
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT3,IERR)
    ENDIF
  ENDIF

!     saturation
  IF( MYPROC == 1 ) NVARID = NCVID(NPOS,'SoilSat',IERR)
  DO JK=1,NCSS
    IF(LEVGEN)THEN
      DO JL=1,NPOI
        ZVALUE(JL)=S3RWSATM3D(JL,JK)
      ENDDO
    ELSE
      ZVALUE(:)=RWSAT
    ENDIF
    CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZVALUE(:),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTCLIM:")
    IF( MYPROC == 1 ) THEN
      ZOUTPUT3(:,JK)= UNPACK(ZBUF,LMASK(ISTP:IENP),RMISS)
    ENDIF
  ENDDO
  IF( MYPROC == 1 ) THEN
    IF(NACCUR == 1)THEN
      ZOUTPUT3S(:,:)=ZOUTPUT3(:,:)
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT3S,IERR)
    ELSE
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT3,IERR)
    ENDIF
  ENDIF

  !     root fraction
  IF( MYPROC == 1 ) NVARID = NCVID(NPOS,'RootDist',IERR)
! !   do jl=0,20
! !     print*,'cov',jl,ZRVCOV(jl)
! !     do jk=1,ncss
! !       print*,'root',jk,ZRVROOTSA(JK,jl)
! !     enddo
! !   enddo
! ! call abort()
! !   print*,ZRVCOV
! !   print*,ZRVROOTSA
! !   print*,'**'
!   do jl=1,npoi
!     ITYPL=NINT(VFTVL(JL))
!     ITYPH=NINT(VFTVH(JL))
!     ZCVL=VFCVL(JL)*ZRVCOV(ITYPL)
!     ZCVH=VFCVH(JL)*ZRVCOV(ITYPH)
!     print*,jl,ITYPL,ZRVCOV(ITYPL),ZCVL
!     print*,jl,ITYPH,ZRVCOV(ITYPH),ZCVH
!   enddo
  DO JK=1,NCSS
!     print*,'jk1 ',jk
    ZVALUE(:) = 0._JPRB
    DO JL=1,NPOI
      ITYPL=NINT(VFTVL(JL))
      ITYPH=NINT(VFTVH(JL))
      ZCVL=VFCVL(JL)*S2RVCOVL2D(JL)
      ZCVH=VFCVH(JL)*S2RVCOVH2D(JL)
      IF((ZCVL > 0._JPRB .AND. S3RVROOTSAL3D(JL,JK) > 0._JPRB) .OR.&
         &(ZCVH > 0._JPRB .AND. S3RVROOTSAH3D(JL,JK) > 0._JPRB)) THEN
        ZVALUE(JL)=(ZCVL*S3RVROOTSAL3D(JL,JK)+ZCVH*S3RVROOTSAH3D(JL,JK))/(ZCVL+ZCVH)
      ENDIF
    ENDDO
!      print*,'jk2 ',jk
    CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZVALUE(:),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTCLIM:")
    IF( MYPROC == 1 ) THEN
      ZOUTPUT3(:,JK)= UNPACK(ZBUF,LMASK(ISTP:IENP),RMISS)
    ENDIF
  ENDDO
  IF( MYPROC == 1 ) THEN
    IF(NACCUR == 1)THEN
      ZOUTPUT3S(:,:)=ZOUTPUT3(:,:)
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT3S,IERR)
    ELSE
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT3,IERR)
    ENDIF
  ENDIF
ENDIf

DEALLOCATE (ZOUTPUT3)  !INFLATED OUTPUT ARRAY
DEALLOCATE (ZOUTPUT3S) !INFLATED OUTPUT ARRAY
DEALLOCATE (ZVALUE)
DEALLOCATE (ZBUF)

IF (LHOOK) CALL DR_HOOK('WRTCLIM',1,ZHOOK_HANDLE)

RETURN
END SUBROUTINE WRTCLIM
